rm(list = ls())
library(foreign)
library(stargazer)
library(survival)
library(Amelia)

set.seed(981237)

#set working directory
source('functions.R')
source('pool_compare.R')
load('main_data.Rdata')

#----------------------------------------------------------------------------------------
#Table A1
#----------------------------------------------------------------------------------------
sumStatVars_spells <- c('cus_sr', 'e_migdppclnL2', 'e_migdpgroL2', 'ln.tpopL2', 'ln.oil_gas_valuePOP_2014L2', 'partB', 'persB', 'mil', 'mon', 'sv_comm', 'InstalledCom')
corSpells <- cor(cus_spells[,sumStatVars_spells], , use = 'complete.obs')
colnames(corSpells) <- c('Rev', 'GDP','Growth', 'Pop', 'Oil/Gas', 'Party', 'Pers.', 'Mil.', 'Mon.', 'Communist', 'Foreign installed communist')
rownames(corSpells) <- c('Rev', 'GDP','Growth', 'Pop', 'Oil/Gas', 'Party', 'Pers.', 'Mil.', 'Mon.', 'Communist', 'Foreign installed communist')

TableA1 <- stargazer(corSpells, title = "Correlation Table, Fixed Covariates", notes.append = FALSE, align = TRUE, float = FALSE)
for (entry in dictSummstats){
  TableA1 <- gsub(entry, names(dictSummstats)[dictSummstats == entry], TableA1)
}

#----------------------------------------------------------------------------------------
#Table A2
#----------------------------------------------------------------------------------------
#Cox Regressions With Time Varying Covariates
form1.tv <- Surv(cus_t, cus_tplus1, cus_fail) ~ cus_sr + e_migdppclnL1 + e_migdpgroL1 + ln.tpopL1
model1.tv.base <- ameliaCoxEstimate(ameliaObject = d_mi_ts, formul = update(form1.tv, . ~ . - cus_sr))
model1.tv <- ameliaCoxEstimate(ameliaObject = d_mi_ts, formul = form1.tv)
model1.tv$summary

form2.tv <- Surv(cus_t, cus_tplus1, cus_fail) ~ cus_sr + e_migdppclnL1 + e_migdpgroL1 + ln.tpopL1 + ln.oil_gas_valuePOP_2014L1
model2.tv.base <- ameliaCoxEstimate(ameliaObject = d_mi_ts, formul = update(form2.tv, . ~ . - cus_sr))
model2.tv <- ameliaCoxEstimate(ameliaObject = d_mi_ts, formul = form2.tv)
model2.tv$summary

form3.tv <- Surv(cus_t, cus_tplus1, cus_fail) ~ cus_sr + e_migdppclnL1 + e_migdpgroL1+ ln.tpopL1 + ln.oil_gas_valuePOP_2014L1  + strata(e_regionpol)
model3.tv.base <- ameliaCoxEstimate(ameliaObject = d_mi_ts, formul = update(form3.tv, . ~ . -cus_sr))
model3.tv <- ameliaCoxEstimate(ameliaObject = d_mi_ts, formul = form3.tv)
model3.tv$summary

form4.tv <- Surv(cus_t, cus_tplus1, cus_fail) ~ cus_sr + e_migdppclnL1 + e_migdpgroL1 + ln.tpopL1 + ln.oil_gas_valuePOP_2014L1 + strata(cus_ARevnum)
model4.tv.base <- ameliaCoxEstimate(ameliaObject = d_mi_ts, formul = update(form4.tv, . ~ . -cus_sr))
model4.tv <- ameliaCoxEstimate(ameliaObject = d_mi_ts, formul = form4.tv)
model4.tv$summary

form5.tv <- Surv(cus_t, cus_tplus1, cus_fail) ~ cus_sr+ e_migdppclnL1 + e_migdpgroL1 + ln.tpopL1 + ln.oil_gas_valuePOP_2014L1 + partB +persB + mil + mon
model5.tv.base <- ameliaCoxEstimate(ameliaObject = d_mi_ts, formul = update(form5.tv, . ~ . -cus_sr))
model5.tv <- ameliaCoxEstimate(ameliaObject = d_mi_ts, formul = form5.tv)
model5.tv$summary

form6.tv <- Surv(cus_t, cus_tplus1, cus_fail) ~ cus_sr + e_migdppclnL1 + e_migdpgroL1 + ln.tpopL1 + ln.oil_gas_valuePOP_2014L1 + partB +persB + mil + mon  + strata(cus_ARevnum)
model6.tv.base <- ameliaCoxEstimate(ameliaObject = d_mi_ts, formul = update(form6.tv, . ~ . -cus_sr))
model6.tv <- ameliaCoxEstimate(ameliaObject = d_mi_ts, formul = form6.tv)
model6.tv$summary

form7.tv <- Surv(cus_t, cus_tplus1, cus_fail) ~ cus_sr + e_migdppclnL1 + e_migdpgroL1 + ln.tpopL1  + ln.oil_gas_valuePOP_2014L1+ partB +persB + mil + mon + sv_comm + strata(cus_ARevnum)
model7.tv.base <- ameliaCoxEstimate(ameliaObject = d_mi_ts, formul = update(form7.tv, . ~ . -cus_sr))
model7.tv <- ameliaCoxEstimate(ameliaObject = d_mi_ts, formul = form7.tv)
model7.tv$summary

form8.tv <- Surv(cus_t, cus_tplus1, cus_fail) ~ cus_sr + e_migdppclnL1 + e_migdpgroL1 + ln.tpopL1 +ln.oil_gas_valuePOP_2014L1  + sv_comm + partB +persB + mil + mon + InstalledCom + strata(cus_ARevnum)
model8.tv.base <- ameliaCoxEstimate(ameliaObject = d_mi_ts, formul = update(form8.tv, . ~ . -cus_sr))
model8.tv <- ameliaCoxEstimate(ameliaObject = d_mi_ts, formul = form8.tv)
model8.tv$summary

ref.models.tv <- list(
  model1.tv.base,
  model2.tv.base,
  model3.tv.base,
  model4.tv.base,
  model5.tv.base,
  model6.tv.base,
  model7.tv.base,
  model8.tv.base
)

modelListObject.tv <-
  CoxmodelList(
    model1.tv,
    model2.tv,
    model3.tv,
    model4.tv,
    model5.tv,
    model6.tv,
    model7.tv,
    model8.tv
  )

ModelLabels <- c('Non-stratified',
                 'Stratified',
                 'Controlling for regime type (post-treatment)',
                 'Controlling for communism (post-treatment)')

column.separate <- c(2, 2, 2, 2)

TableA2 <- CoxTable(modelListObject = modelListObject.tv,
                      ref.Models = ref.models.tv,
                      dict = dictSummstatsL1,
                      column.labels = ModelLabels,
                      column.separate = column.separate,
                      dep.var.labels.include = FALSE,
                      align = TRUE,
                      float = FALSE,
                      font.size = 'footnotesize')


#----------------------------------------------------------------------------------------
#Table A3
#----------------------------------------------------------------------------------------
#Listwise deletion, time varying covariates
dat$cus_tplus1 <- dat$cus_t + 1
dset <- subset(dat, cus_AR%in%1 & cus_syear >= 1900)
dset$id <- factor(1:nrow(dset))

form1.lw <- Surv(cus_t, cus_tplus1, cus_fail) ~ cus_sr + e_migdppclnL1 + e_migdpgroL1  + ln.tpopL1 + cluster(id)
model1.lw.base <- coxph(update(form1.lw, . ~ . - cus_sr - cluster(id)), data = dset)
model1.lw <- coxph(form1.lw, data = dset)
model1.lw.anov <- coxph(update(form1.lw, . ~ . -cluster(id)), data = dset)
model1.lw.lm <- lm(update(form1.lw, cus_fail ~ . - cluster(id)), data = dset)
#summary(model1.tv)

form2.lw <- Surv(cus_t, cus_tplus1, cus_fail) ~ cus_sr + e_migdppclnL1 + e_migdpgroL1  + ln.tpopL1  + ln.oil_gas_valuePOP_2014L1 + cluster(id)
model2.lw.base <- coxph(update(form2.lw, . ~ . - cus_sr - cluster(id)), data = dset)
model2.lw <- coxph(form2.lw, data = dset)
model2.lw.anov <- coxph(update(form2.lw, . ~ . -cluster(id)), data = dset)
model2.lw.lm <- lm(update(form2.lw, cus_fail ~ . - cluster(id)), data = dset)
#summary(form2.lw)

form3.lw <- Surv(cus_t, cus_tplus1, cus_fail) ~ cus_sr + e_migdppclnL1 + e_migdpgroL1  + ln.tpopL1  + ln.oil_gas_valuePOP_2014L1 + strata(e_regionpol) + cluster(id)
model3.lw.base <- coxph(update(form3.lw, . ~ . - cus_sr - cluster(id)), data = dset)
model3.lw <- coxph(form3.lw, data = dset)
model3.lw.anov <- coxph(update(form3.lw, . ~ . -cluster(id)), data = dset)
model3.lw.lm <- lm(update(form3.lw, cus_fail ~ . - cluster(id)), data = dset)
#summary(model3.lw)

form4.lw <- Surv(cus_t, cus_tplus1, cus_fail) ~ cus_sr + e_migdppclnL1 + e_migdpgroL1  + ln.tpopL1  + ln.oil_gas_valuePOP_2014L1 + strata(cus_ARevnum) + cluster(id)
model4.lw.base <- coxph(update(form4.lw, . ~ . - cus_sr - cluster(id)), data = dset)
model4.lw <- coxph(form4.lw, data = dset)
model4.lw.anov <- coxph(update(form4.lw, . ~ . -cluster(id)), data = dset)
model4.lw.lm <- lm(update(form4.lw, cus_fail ~ . - cluster(id)), data = dset)
summary(model4.lw)

form5.lw <- Surv(cus_t, cus_tplus1, cus_fail) ~ cus_sr + e_migdppclnL1 + e_migdpgroL1  + ln.tpopL1  + ln.oil_gas_valuePOP_2014L1 + strata(e_regionpol) + partB + persB + mil + mon + cluster(id)
model5.lw.base <- coxph(update(form5.lw, . ~ . - cus_sr - cluster(id)), data = dset)
model5.lw <- coxph(form5.lw, data = dset)
model5.lw.anov <- coxph(update(form5.lw, . ~ . -cluster(id)), data = dset)
model5.lw.lm <- lm(update(form5.lw, cus_fail ~ . - cluster(id)), data = dset)
summary(model5.lw)

form6.lw <- Surv(cus_t, cus_tplus1, cus_fail) ~ cus_sr + e_migdppclnL1 + e_migdpgroL1  + ln.tpopL1  + ln.oil_gas_valuePOP_2014L1 + strata(cus_ARevnum) + partB + persB + mil + mon + cluster(id)
model6.lw.base <- coxph(update(form6.lw, . ~ . - cus_sr - cluster(id)), data = dset)
model6.lw <- coxph(form6.lw, data = dset)
model6.lw.anov <- coxph(update(form6.lw, . ~ . -cluster(id)), data = dset)
model6.lw.lm <- lm(update(form6.lw, cus_fail ~ . - cluster(id)), data = dset)
summary(model6.lw)

form7.lw <- Surv(cus_t, cus_tplus1, cus_fail) ~ cus_sr + e_migdppclnL1 + e_migdpgroL1  + ln.tpopL1  + ln.oil_gas_valuePOP_2014L1 + strata(cus_ARevnum) + partB + persB + mil + mon + sv_comm + cluster(id)
model7.lw.base <- coxph(update(form7.lw, . ~ . - cus_sr - cluster(id)), data = dset)
model7.lw <- coxph(form7.lw, data = dset)
model7.lw.anov <- coxph(update(form7.lw, . ~ . -cluster(id)), data = dset)
model7.lw.lm <- lm(update(form7.lw, cus_fail ~ . - cluster(id)), data = dset)
summary(model7.lw)

form8.lw <- Surv(cus_t, cus_tplus1, cus_fail) ~ cus_sr + e_migdppclnL1 + e_migdpgroL1  + ln.tpopL1  + ln.oil_gas_valuePOP_2014L1 + strata(cus_ARevnum) + partB + persB + mil + mon + sv_comm + InstalledCom + cluster(id)
model8.lw.base <- coxph(update(form8.lw, . ~ . - cus_sr - cluster(id)), data = dset)
model8.lw <- coxph(form8.lw, data = dset)
model8.lw.anov <- coxph(update(form8.lw, . ~ . -cluster(id)), data = dset)
model8.lw.lm <- lm(update(form8.lw, cus_fail ~ . - cluster(id)), data = dset)
summary(model8.lw)

coefs.lw <- list(
  model1.lw$coefficients,
  model2.lw$coefficients,
  model3.lw$coefficients,
  model4.lw$coefficients,
  model5.lw$coefficients,
  model6.lw$coefficients,
  model7.lw$coefficients,
  model8.lw$coefficients
)

ses.lw <- list(
  summary(model1.lw)$coefficients[, 4],
  summary(model2.lw)$coefficients[, 4],
  summary(model3.lw)$coefficients[, 4],
  summary(model4.lw)$coefficients[, 4],
  summary(model5.lw)$coefficients[, 4],
  summary(model6.lw)$coefficients[, 4],
  summary(model7.lw)$coefficients[, 4],
  summary(model8.lw)$coefficients[, 4]
)

nbreakdowns <- c(
  sum(model1.lw$y[,3]),
  sum(model2.lw$y[,3]),
  sum(model3.lw$y[,3]),
  sum(model4.lw$y[,3]),
  sum(model5.lw$y[,3]),
  sum(model6.lw$y[,3]),
  sum(model7.lw$y[,3]),
  sum(model8.lw$y[,3])
)
nobs <- c(
  length(model1.lw$y[,3]),
  length(model2.lw$y[,3]),
  length(model3.lw$y[,3]),
  length(model4.lw$y[,3]),
  length(model5.lw$y[,3]),
  length(model6.lw$y[,3]),
  length(model7.lw$y[,3]),
  length(model8.lw$y[,3])
  
)

logliktests <- c(
  deci(anova(model1.lw.anov, model1.lw.base)$'P(>|Chi|)'[2], 2),
  deci(anova(model2.lw.anov, model2.lw.base)$'P(>|Chi|)'[2], 2),
  deci(anova(model3.lw.anov, model2.lw.base)$'P(>|Chi|)'[2], 2),
  deci(anova(model4.lw.anov, model4.lw.base)$'P(>|Chi|)'[2], 2),
  deci(anova(model5.lw.anov, model5.lw.base)$'P(>|Chi|)'[2], 2),
  deci(anova(model6.lw.anov, model5.lw.base)$'P(>|Chi|)'[2], 2),
  deci(anova(model7.lw.anov, model5.lw.base)$'P(>|Chi|)'[2], 2),
  deci(anova(model8.lw.anov, model5.lw.base)$'P(>|Chi|)'[2], 2)
)

ModelLabels <- c('Non-stratified',
                 'Stratified',
                 'Controlling for regime type (post-treatment)',
                 'Controlling for communism (post-treatment)'
)

column.separate <- c(2, 2, 2, 2)
TableA3 <- stargazer(
  model1.lw.lm,
  model2.lw.lm,
  model3.lw.lm,
  model4.lw.lm,
  model5.lw.lm,
  model6.lw.lm,
  model7.lw.lm,
  model8.lw.lm,
  coef = coefs.lw,
  se = ses.lw,
  column.separate = column.separate,
  column.labels = ModelLabels,
  dep.var.labels = c('', '', '', '', '', '', '', ''),
  dep.var.caption= '',
  align = TRUE,
  float = FALSE,
  font.size = 'footnotesize',
  omit = c('strata', 'Constant', 'cluster'),
  omit.stat = c('adj.rsq', 'ser', 'f', 'rsq', 'll', 'n'),
  add.lines = list(
                    c('Stratified by geographic region', 'No', 'No', "Yes", 'No', "No", "No", "No", "No"),
                    c('Stratified by event number', 'No', 'No', 'No', "Yes", "Yes", "Yes", "Yes", "Yes"),
                    c('N observations', nobs), 
                    c('N breakdowns', nbreakdowns),
                    c('Multiple imputation', 'No', 'No', 'No', 'No', 'No', 'No', 'No', 'No'),
                    c('$p$-value of significance test:', logliktests) )
)

for (entry in dictSummstatsL1){
  TableA3 <- gsub(entry, names(dictSummstatsL1)[dictSummstatsL1 == entry], TableA3)
}


#----------------------------------------------------------------------------------------
#Table A4
#----------------------------------------------------------------------------------------
#Discard long-lasting rev regimes
#Table A4 a)
D_1 <- D_2 <- D_3 <- D_4 <- D_5 <- D_6 <- d_mi_spells_cus
D_1$imputations <- lapply(d_mi_spells_cus$imputations,
                          function(x) subset(x, !(cus_caseid %in% 'Mexico 1915-2000')))
D_2$imputations <- lapply(d_mi_spells_cus$imputations,
                          function(x) subset(x, !(cus_caseid %in% c('Mexico 1915-2000', 'Russia 1917-1991'))))
D_3$imputations <- lapply(d_mi_spells_cus$imputations,
                          function(x) subset(x, !(cus_caseid %in% c('Mexico 1915-2000', 'Russia 1917-1991', 'China 1949-'))))
D_4$imputations <- lapply(d_mi_spells_cus$imputations,
                          function(x) return(subset(x, !(cus_caseid %in% c('Mexico 1915-2000', 'Russia 1917-1991', 'China 1949-', 'Vietnam 1954-')))))
D_5$imputations <- lapply(d_mi_spells_cus$imputations,
                          function(x) return(subset(x, !(cus_caseid %in% c('Mexico 1915-2000', 'Russia 1917-1991', 'China 1949-', 'Vietnam 1954-', 'Cuba 1959-')))))
D_6$imputations <- lapply(d_mi_spells_cus$imputations,
                          function(x) return(subset(x, !(cus_caseid %in% c('Mexico 1915-2000', 'Russia 1917-1991', 'China 1949-', 'Vietnam 1954-', 'Cuba 1959-', 'Algeria 1962-')))))

form1 <- Surv(cus_t_surv, cus_fail) ~ cus_sr + e_migdppclnL2 + e_migdpgroL2 + ln.tpopL2 + ln.oil_gas_valuePOP_2014L2
model1.base <- ameliaCoxEstimate(ameliaObject = D_1, formul = update(form1, . ~ . -cus_sr))
model1 <- ameliaCoxEstimate(ameliaObject = D_1, formul = form1)
model1$summary

model2.base <- ameliaCoxEstimate(ameliaObject = D_2, formul = update(form1, . ~ . -cus_sr))
model2 <- ameliaCoxEstimate(ameliaObject = D_2, formul = form1)
model2$summary

model3.base <- ameliaCoxEstimate(ameliaObject = D_3, formul = update(form1, . ~ . -cus_sr))
model3 <- ameliaCoxEstimate(ameliaObject = D_3, formul = form1)
model3$summary

model4.base <- ameliaCoxEstimate(ameliaObject = D_4, formul = update(form1, . ~ . -cus_sr))
model4 <- ameliaCoxEstimate(ameliaObject = D_4, formul = form1)
model4$summary

model5.base <- ameliaCoxEstimate(ameliaObject = D_5, formul = update(form1, . ~ . -cus_sr))
model5 <- ameliaCoxEstimate(ameliaObject = D_5, formul = form1)
model5$summary

model6.base <- ameliaCoxEstimate(ameliaObject = D_6, formul = update(form1, . ~ . -cus_sr))
model6 <- ameliaCoxEstimate(ameliaObject = D_6, formul = form1)
model6$summary

#fit1 <- coxph(form3, data = cus_spells)
#fit0 <- coxph(update(form3, . ~ . - cus_sr), data = cus_spells)
ref.models <- list(
  model1.base,
  model2.base,
  model3.base,
  model4.base,
  model5.base,
  model6.base
)

modelListObject <-
  CoxmodelList(
    model1,
    model2,
    model3,
    model4,
    model5,
    model6
  )
ModelLabels <- c('N$-$1',
                 'N$-$2',
                 'N$-$3',
                 'N$-$4',
                 'N$-$5',
                 'N$-$6'
)

strat.custom <- list(c('Stratified by geographic region', 'No', "No", 'No', 'No',  'No', "No", 'No', "No"),
                     c('Stratified by event number', 'No', 'No', 'No', "No",  "No", "No", "No", "No"))

TableA4a_t <- CoxTable(modelListObject = modelListObject,
                              ref.Models = ref.models,
                              dict = dictSummstats,
                              column.labels = ModelLabels,
                              dep.var.labels = 'Authoritarian Regime Breakdown',
                              align = TRUE,
                              float = FALSE,
                              strat.custom = strat.custom,
                              font.size = 'footnotesize')

TableA4a <- c(TableA4a_t[1:(length(TableA4a_t)-5)], "\\end{tabular} " , "\\endgroup ")


#Table A4 b)
form1 <- Surv(cus_t_surv, cus_fail) ~ cus_sr + e_migdppclnL2 + e_migdpgroL2 + ln.tpopL2 + ln.oil_gas_valuePOP_2014L2 + partB + mil + mon + persB
model1.base <- ameliaCoxEstimate(ameliaObject = D_1, formul = update(form1, . ~ . -cus_sr))
model1 <- ameliaCoxEstimate(ameliaObject = D_1, formul = form1)
model1$summary

model2.base <- ameliaCoxEstimate(ameliaObject = D_2, formul = update(form1, . ~ . -cus_sr))
model2 <- ameliaCoxEstimate(ameliaObject = D_2, formul = form1)
model2$summary

model3.base <- ameliaCoxEstimate(ameliaObject = D_3, formul = update(form1, . ~ . -cus_sr))
model3 <- ameliaCoxEstimate(ameliaObject = D_3, formul = form1)
model3$summary

model4.base <- ameliaCoxEstimate(ameliaObject = D_4, formul = update(form1, . ~ . -cus_sr))
model4 <- ameliaCoxEstimate(ameliaObject = D_4, formul = form1)
model4$summary

model5.base <- ameliaCoxEstimate(ameliaObject = D_5, formul = update(form1, . ~ . -cus_sr))
model5 <- ameliaCoxEstimate(ameliaObject = D_5, formul = form1)
model5$summary

model6.base <- ameliaCoxEstimate(ameliaObject = D_6, formul = update(form1, . ~ . -cus_sr))
model6 <- ameliaCoxEstimate(ameliaObject = D_6, formul = form1)
model6$summary

ref.models <- list(
  model1.base,
  model2.base,
  model3.base,
  model4.base,
  model5.base,
  model6.base
)

modelListObject <-
  CoxmodelList(
    model1,
    model2,
    model3,
    model4,
    model5,
    model6
  )
ModelLabels <- c('N$-$1',
                 'N$-$2',
                 'N$-$3',
                 'N$-$4',
                 'N$-$5',
                 'N$-$6'
)

strat.custom <- list(c('Stratified by geographic region', 'No', "No", 'No', 'No',  'No', "No", 'No', "No"),
                     c('Stratified by event number', 'No', 'No', 'No', "No",  "No", "No", "No", "No"))

TableA4b <- CoxTable(modelListObject = modelListObject,
                                    ref.Models = ref.models,
                                    dict = dictSummstats,
                                    column.labels = ModelLabels,
                                    dep.var.labels = 'Authoritarian Regime Breakdown',
                                    align = TRUE,
                                    float = FALSE,
                                    strat.custom = strat.custom,
                                    font.size = 'footnotesize')




#----------------------------------------------------------------------------------------
#Table A5
#----------------------------------------------------------------------------------------
#Include Finland and Hungary
form1_N20 <- Surv(cus_t_surv, cus_fail) ~ cus_sr + e_migdppclnL2 + e_migdpgroL2 + ln.tpopL2
model1_N20.base <- ameliaCoxEstimate(ameliaObject = d_mi_spells_N20, formul = update(form1_N20, . ~ . - cus_sr))
model1_N20 <- ameliaCoxEstimate(ameliaObject = d_mi_spells_N20, formul = form1_N20)
model1_N20$summary;model1_N20$N;sum(d_mi_spells_N20$imputations[[1]]$cus_sr)
model1_N20$summary; model1_N20$N

form2_N20 <- Surv(cus_t_surv, cus_fail) ~ cus_sr + e_migdppclnL2 + e_migdpgroL2 + ln.tpopL2 + ln.oil_gas_valuePOP_2014L2
model2_N20.base <- ameliaCoxEstimate(ameliaObject = d_mi_spells_N20, formul = update(form2_N20, . ~ . - cus_sr))
model2_N20 <- ameliaCoxEstimate(ameliaObject = d_mi_spells_N20, formul = form2_N20)
model2_N20$summary; model2_N20$N

form3_N20 <- Surv(cus_t_surv, cus_fail) ~ cus_sr + e_migdppclnL2 + e_migdpgroL2+ ln.tpopL2 + ln.oil_gas_valuePOP_2014L2  + strata(e_regionpol)
model3_N20.base <- ameliaCoxEstimate(ameliaObject = d_mi_spells_N20, formul = update(form3_N20, . ~ . -cus_sr))
model3_N20 <- ameliaCoxEstimate(ameliaObject = d_mi_spells_N20, formul = form3_N20, DV_Label = 'covFixed')
model3_N20$summary; model3_N20$N

form4_N20 <- Surv(cus_t_surv, cus_fail) ~ cus_sr + e_migdppclnL2 + e_migdpgroL2 + ln.tpopL2 + ln.oil_gas_valuePOP_2014L2 + strata(cus_ARevnum)
model4_N20.base <- ameliaCoxEstimate(ameliaObject = d_mi_spells_N20, formul = update(form4_N20, . ~ . -cus_sr))
model4_N20 <- ameliaCoxEstimate(ameliaObject = d_mi_spells_N20, formul = form4_N20)
model4_N20$summary; model4_N20$N


form5_N20 <- Surv(cus_t_surv, cus_fail) ~ cus_sr+ e_migdppclnL2 + e_migdpgroL2 + ln.tpopL2 + ln.oil_gas_valuePOP_2014L2 + partB +persB + mil + mon
model5_N20.base <- ameliaCoxEstimate(ameliaObject = d_mi_spells_N20, formul = update(form5_N20, . ~ . -cus_sr))
model5_N20 <- ameliaCoxEstimate(ameliaObject = d_mi_spells_N20, formul = form5_N20)
model5_N20$summary; model5_N20$N

form6_N20 <- Surv(cus_t_surv, cus_fail) ~ cus_sr + e_migdppclnL2 + e_migdpgroL2 + ln.tpopL2 + ln.oil_gas_valuePOP_2014L2 + partB +persB + mil + mon  + strata(cus_ARevnum)
model6_N20.base <- ameliaCoxEstimate(ameliaObject = d_mi_spells_N20, formul = update(form6_N20, . ~ . -cus_sr))
model6_N20 <- ameliaCoxEstimate(ameliaObject = d_mi_spells_N20, formul = form6_N20)
model6_N20$summary; model6_N20$N

form7_N20 <- Surv(cus_t_surv, cus_fail) ~ cus_sr + e_migdppclnL2 + e_migdpgroL2 + ln.tpopL2  + ln.oil_gas_valuePOP_2014L2+ partB +persB + mil + mon + sv_comm+ strata(cus_ARevnum)
model7_N20.base <- ameliaCoxEstimate(ameliaObject = d_mi_spells_N20, formul = update(form7_N20, . ~ . -cus_sr))
model7_N20 <- ameliaCoxEstimate(ameliaObject = d_mi_spells_N20, formul = form7_N20)
model7_N20$summary; model7_N20$N

form8_N20 <- Surv(cus_t_surv, cus_fail) ~ cus_sr + e_migdppclnL2 + e_migdpgroL2 + ln.tpopL2 +ln.oil_gas_valuePOP_2014L2  + sv_comm+ partB +persB + mil + mon + InstalledCom+ strata(cus_ARevnum)
model8_N20.base <- ameliaCoxEstimate(ameliaObject = d_mi_spells_N20, formul = update(form8_N20, . ~ . -cus_sr))
model8_N20 <- ameliaCoxEstimate(ameliaObject = d_mi_spells_N20, formul = form8_N20)
model8_N20$summary;model8_N20$N

ref.models.N20 <- list(
  model1_N20.base,
  model2_N20.base,
  model3_N20.base,
  model4_N20.base,
  model5_N20.base,
  model6_N20.base,
  model7_N20.base,
  model8_N20.base
)

modelListObject.N20 <-
  CoxmodelList(
    model1_N20,
    model2_N20,
    model3_N20,
    model4_N20,
    model5_N20,
    model6_N20,
    model7_N20,
    model8_N20
  )

ModelLabels <- c('Non-stratified',
                 'Stratified',
                 'Controlling for regime type (post-treatment)',
                 'Controlling for communism (post-treatment)')

column.separate <- c(2, 2, 2, 2)

TableA5 <- CoxTable(modelListObject = modelListObject.N20,
                        ref.Models = ref.models.N20,
                        dict = dictSummstatsL2,
                        column.labels = ModelLabels,
                        column.separate = column.separate,
                        dep.var.labels.include = FALSE,
                        align = TRUE,
                        float = FALSE,
                        font.size = 'footnotesize')


#----------------------------------------------------------------------------------------
#Table A6
#----------------------------------------------------------------------------------------
#Bootstrap
library(boot)
set.seed(4132478)
print('bootstrap1')
form1_btsp <- Surv(cus_t_surv, cus_fail) ~ cus_sr + partB +persB + mil + mon 
model1_btsp.base <- coxph(update(form1_btsp, . ~ . - cus_sr), data = cus_spells)
model1_btsp <- coxph(form1_btsp, data = cus_spells)
bootstrap1.out <- boot(cus_spells, coxph_boot, R = 100000, frmla = form1_btsp)
btsp_se1 <- apply(bootstrap1.out$t, 2, sd); names(btsp_se1) <- names(model1_btsp$coefficients)

print('bootstrap2')
form2_btsp <- Surv(cus_t_surv, cus_fail) ~ cus_sr + partB +persB + mil + mon+ sv_comm
model2_btsp.base <- coxph(update(form2_btsp, . ~ . - cus_sr), data = cus_spells)
model2_btsp <- coxph(form2_btsp, data = cus_spells)
bootstrap2.out <- boot(cus_spells, coxph_boot, R = 100000, frmla = form2_btsp)
btsp_se2 <- apply(bootstrap2.out$t, 2, sd); names(btsp_se2) <- names(model2_btsp$coefficients)
summary(model2_btsp)

print('bootstrap3')
form3_btsp <- Surv(cus_t_surv, cus_fail) ~ cus_sr + partB +persB + mil + mon+ sv_comm + InstalledCom
model3_btsp.base <- coxph(update(form2_btsp, . ~ . - cus_sr), data = cus_spells)
model3_btsp <- coxph(form3_btsp, data = cus_spells)
bootstrap3.out <- boot(cus_spells, coxph_boot, R = 10000, frmla = form3_btsp)
btsp_se3 <- apply(bootstrap3.out$t, 2, function(x) sd(x, na.rm=T)); names(btsp_se3) <- names(model3_btsp$coefficients)

coefs.btsp <- list(
  model1_btsp$coefficients,
  model2_btsp$coefficients,
  model3_btsp$coefficients
)

ses.btsp <- list(
  btsp_se1,
  btsp_se2,
  btsp_se3
)

nbreakdowns <- c(
  sum(model1_btsp$y[,2]),
  sum(model2_btsp$y[,2]),
  sum(model3_btsp$y[,2])
)

nobs <- c(
  length(model1_btsp$y[,2]),
  length(model2_btsp$y[,2]),
  length(model3_btsp$y[,2])
  )

logliktests <- c(
  deci(anova(model1_btsp, model1_btsp.base)$'P(>|Chi|)'[2], 3),
  deci(anova(model2_btsp, model2_btsp.base)$'P(>|Chi|)'[2], 3),
  deci(anova(model3_btsp, model3_btsp.base)$'P(>|Chi|)'[2], 3)
)

column.separate <- c(1)
TableA6 <- stargazer(
  model1_btsp,
  model2_btsp,
  model3_btsp,
  coef = coefs.btsp,
  se = ses.btsp,
  column.separate = column.separate,
  dep.var.labels = c(''),
  dep.var.caption= '',
  align = TRUE,
  float = FALSE,
  font.size = 'footnotesize',
  omit = c('strata', 'Constant', 'cluster'),
  omit.stat = c('all'),
  add.lines = list(c('N observations', nobs), c('N breakdowns', nbreakdowns),
                   c('Stratified by geographic region', 'No', 'No', "No"),
                   c('Stratified by event number', 'No', 'No', 'No'),
                   c('$p$-value of significance test:', logliktests) )
)
for (entry in dictSummstatsL1){
  TableA6 <- gsub(entry, names(dictSummstatsL1)[dictSummstatsL1 == entry], TableA6)
}


#----------------------------------------------------------------------------------------
#Table A7
#----------------------------------------------------------------------------------------
#Stratify by country

form1 <- Surv(cus_t_surv, cus_fail) ~ cus_sr + e_migdppclnL2 + e_migdpgroL2 + ln.tpopL2 +  strata(ccode)
model1.base <- ameliaCoxEstimate(ameliaObject = d_mi_spells_cus, formul = update(form1, . ~ . -cus_sr))
model1 <- ameliaCoxEstimate(ameliaObject = d_mi_spells_cus, formul = form1)
model1$summary

form2 <- Surv(cus_t_surv, cus_fail) ~ cus_sr + e_migdppclnL2 + e_migdpgroL2 + ln.tpopL2 + ln.oil_gas_valuePOP_2014L2  + strata(ccode) 
model2.base <- ameliaCoxEstimate(ameliaObject = d_mi_spells_cus, formul = update(form2, . ~ . -cus_sr))
model2 <- ameliaCoxEstimate(ameliaObject = d_mi_spells_cus, formul = form2)
model2$summary

ref.models <- list(
  model1.base,
  model2.base
)

modelListObject <-
  CoxmodelList(
    model1,
    model2
  )

strat.custom <- list(c('Stratified by geographic region', 'No', "No"),
                     c('Stratified by event number', 'No', 'No'))

TableA7 <- CoxTable(
  modelListObject = modelListObject,
  ref.Models = ref.models,
  dict = dictSummstats,
  dep.var.labels.include = FALSE,
  dep.var.labels = 'Authoritarian Regime Breakdown',
  column.labels = '',
  align = TRUE,
  float = FALSE,
  font.size = 'footnotesize',
  strat.custom = strat.custom)


#----------------------------------------------------------------------------------------
#Table A8
#----------------------------------------------------------------------------------------
#Cox regression of coups
form_coups <- Surv(coup_t, coup_tplus1, coup) ~ cus_sr + e_migdppclnL1 + e_migdpgroL1 + ln.tpopL1 + ln.oil_gas_valuePOP_2014L1
model_coup.base <- ameliaCoxEstimate(ameliaObject = d_mi_ts_coups, formul = update(form_coups, . ~ . -cus_sr ))
model_coup <- ameliaCoxEstimate(ameliaObject = d_mi_ts_coups, formul = form_coups)
model_coup$summary

ref.models <- list(
  model_coup.base
)

modelListObject <-
  CoxmodelList(
    model_coup
  )

#column.separate <- c(1, 1)
TableA8_t <- CoxTable(modelListObject = modelListObject,
                        ref.Models = ref.models,
                        dict = dictSummstats,
                        align = TRUE,
                        float = FALSE,
                        font.size = 'footnotesize')

TableA8 <- gsub('breakdowns', 'events', TableA8_t)


#----------------------------------------------------------------------------------------
#Figures A1-A3
#----------------------------------------------------------------------------------------

spells.imp.in <- d_mi_spells_cus
recode_NK <- function(spells.imp.in){
  spells.imp.out <- spells.imp.in
  spells.imp.out$imputations <- lapply(spells.imp.out$imputations, function(dd){
    dd[dd$cus_caseid=='Korea, North 1948-', 'cus_sr'] <- 1
    return(dd)
  })
  return(spells.imp.out)
}
recode_SY <- function(spells.imp.in){
  spells.imp.out <- spells.imp.in
  spells.imp.out$imputations <- lapply(spells.imp.out$imputations, function(dd){
    dd[dd$cus_caseid=='Yemen, South 1967-1990', 'cus_sr'] <- 1
    return(dd)
  })
  return(spells.imp.out)
}
recode_RW <- function(spells.imp.in){
  spells.imp.out <- spells.imp.in
  spells.imp.out$imputations <- lapply(spells.imp.out$imputations, function(dd){
    dd[dd$cus_caseid=='Rwanda 1994-', 'cus_sr'] <- 0
    return(dd)
  })
  return(spells.imp.out)
}
recode_ZM <- function(spells.imp.in){
  spells.imp.out <- spells.imp.in
  spells.imp.out$imputations <- lapply(spells.imp.out$imputations, function(dd){
    dd[dd$cus_caseid=='Zimbabwe 1980-', 'cus_sr'] <- 1
    return(dd)
  })
  return(spells.imp.out)
}

dd <- cus_spells
recode_AL <- function(spells.imp.in){
  spells.imp.out <- spells.imp.in
  spells.imp.out$imputations <- lapply(spells.imp.out$imputations, function(dd){
    newrow <-  dd[dd$cus_caseid%in%'Algeria 1962-', ]
    dd[dd$cus_caseid%in%'Algeria 1962-', 'cus_t_surv'] <- 30
    dd[dd$cus_caseid%in%'Algeria 1962-', 'cus_fail'] <- 1
    dd[dd$cus_caseid%in%'Algeria 1962-', 'cus_caseid'] <- 'Algeria 1962-1992'
    newrow$cus_caseid <- 'Algeria 1992-'
    newrow$cus_t_surv <- 23
    newrow$cus_fail <- 0 
    newrow$e_migdppclnL2 <- dat[dat$Year%in%1991&dat$Country == 'Algeria', 'e_migdppcln']
    newrow$e_migdpgroL2 <- dat[dat$Year%in%1991&dat$Country == 'Algeria', 'e_migdpgro']
    newrow$ln.tpopL2 <- log(dat[dat$Year%in%1991&dat$Country == 'Algeria', 'tpop']+1)
    newrow$ln.oil_gas_valuePOP_2014L2 <- log(dat[dat$Year%in%1991 & dat$Country %in% 'Algeria', 'oil_gas_valuePOP_2014']+1)
    newrow$cus_sr <- 0
    newrow$mil <- 1
    newrow$partB <-0
    newrow$cus_ARevnum <- newrow$cus_ARevnum+1
    dd <- rbind(dd, newrow)
    return(dd)
  })
  return(spells.imp.out)
}

recode_dataset <- function(HU.FI, NK, SY, RW, ZM, AL){
  if (HU.FI == 0){
    dat.out <- d_mi_spells_cus
  }else{
    dat.out <- d_mi_spells_N20
  }
  if (NK == 1){
    dat.out <- recode_NK(dat.out)
  }
  if (SY == 1){
    dat.out <- recode_SY(dat.out)
  }
  if (RW == 1){
    dat.out <- recode_RW(dat.out)
  }
  if (ZM == 1){
    dat.out <- recode_ZM(dat.out)
  }
  if (AL == 1){
    dat.out <- recode_AL(dat.out)
  }
  return(dat.out)
}

form2 <- Surv(cus_t_surv, cus_fail) ~ cus_sr + e_migdppclnL2 + e_migdpgroL2 + ln.tpopL2 + ln.oil_gas_valuePOP_2014L2
form6 <- Surv(cus_t_surv, cus_fail) ~ cus_sr + e_migdppclnL2 + e_migdpgroL2 + ln.tpopL2 + ln.oil_gas_valuePOP_2014L2 + partB+persB+mil+mon + strata(cus_ARevnum)
form8 <- Surv(cus_t_surv, cus_fail) ~ cus_sr + e_migdppclnL2 + e_migdpgroL2 + ln.tpopL2 + ln.oil_gas_valuePOP_2014L2 + partB+persB+mil+mon +sv_comm +InstalledCom + strata(cus_ARevnum)

df <- data.frame(name = rep(NA, 2^6-1),
                 Est = rep(NA, 2^6-1)
)

predeffs <- list(df, df, df)

i <- 1
mod.lab <- c()
for(HU.FI in c(0,1)){
  for(NK in c(0,1)){
    for(SY in c(0,1)){
      for(ZM in c(0,1)){
        for(RW in c(0,1)){
          for(AL in c(0,1)){
              print(c(HU.FI, NK, SY, RW, ZM, AL))
              if (sum(c(HU.FI, NK, SY, RW, ZM, AL) > 0)){
              thisd <- recode_dataset(HU.FI, NK, SY, RW, ZM, AL)
              thismod <- ameliaCoxEstimate(ameliaObject = thisd, formul = form2)
              predeffs[[1]]$Est[i] <- thismod$summary$beta[1]
              predeffs[[1]]$SD[i] <- thismod$summary$se[1]
              thismod <- ameliaCoxEstimate(ameliaObject = thisd, formul = form6)
              predeffs[[2]]$Est[i] <- thismod$summary$beta[1]
              predeffs[[2]]$SD[i] <- thismod$summary$se[1]
              thismod <- ameliaCoxEstimate(ameliaObject = thisd, formul = form8)
              predeffs[[3]]$Est[i] <- thismod$summary$beta[1]
              predeffs[[3]]$SD[i] <- thismod$summary$se[1]
              mod.lab <- c()
              if (NK == 1){mod.lab <- c(mod.lab, 'NK')}
              if (SY == 1){mod.lab <- c(mod.lab, 'SY')}
              if (RW == 1){mod.lab <- c(mod.lab, 'RW')}
              if (ZM == 1){mod.lab <- c(mod.lab, 'ZM')}
              if (AL == 1){mod.lab <- c(mod.lab, 'AL')}
              if (HU.FI == 1){mod.lab <- c(mod.lab, 'HU/FI')}
              print(paste(mod.lab, collapse = '+'))
              predeffs[[1]]$name[i] <- paste(mod.lab, collapse = '+')
              predeffs[[2]]$name[i] <- paste(mod.lab, collapse = '+')
              predeffs[[3]]$name[i] <- paste(mod.lab, collapse = '+')
              i <- i+1
              }
            }
          }
        }
      }
    }
  }


recode.dat <- lapply(predeffs, function(y){
  y$low95 <- y$Est - 1.96*y$SD
  y$low90 <- y$Est - qnorm(0.95)*y$SD
  y$upp90 <- y$Est + qnorm(0.95)*y$SD
  y$upp95 <- y$Est + 1.96*y$SD
  y$name <- factor(as.character(y$name))
  return(y)
})

library(ggplot2)
a1 <- ggplot2::ggplot(recode.dat[[1]], aes(x = Est, y = name)) +
  geom_point(cex = 2) + xlim(-2, .5) +
  geom_errorbarh(aes(xmin = low90, xmax = upp90), height = 0, size = .99, color = 'black') +
  geom_errorbarh(aes(xmin = low95, xmax = upp95), height = 0, size = .39, color = 'black') +
  xlab("") + theme(legend.position="none") +
  ylab("Effect of Revolution on Risk of Breakdown when using alternative codings") +
  geom_vline(xintercept = 0)
ggsave("FigureA1.pdf", a1, width = 6, height = 10)

a2 <- ggplot2::ggplot(recode.dat[[2]], aes(x = Est, y = name)) +
  geom_point(cex = 2) + xlim(-2, .5) +
  geom_errorbarh(aes(xmin = low90, xmax = upp90), height = 0, size = .99, color = 'black') +
  geom_errorbarh(aes(xmin = low95, xmax = upp95), height = 0, size = .39, color = 'black') +
  xlab("") + theme(legend.position="none") +
  ylab("Effect of Revolution on Risk of Breakdown when using alternative codings") +
  geom_vline(xintercept = 0)
ggsave("FigureA2.pdf", a2, width = 6, height = 10)

a3 <- ggplot2::ggplot(recode.dat[[3]], aes(x = Est, y = name)) +
  geom_point(cex = 2) + xlim(-2, .5) +
  geom_errorbarh(aes(xmin = low90, xmax = upp90), height = 0, size = .99, color = 'black') +
  geom_errorbarh(aes(xmin = low95, xmax = upp95), height = 0, size = .39, color = 'black') +
  xlab("") + theme(legend.position="none") +
  ylab("Effect of Revolution on Risk of Breakdown when using alternative codings") +
  geom_vline(xintercept = 0)
ggsave("FigureA3.pdf", a3, width = 6, height = 10)


#----------------------------------------------------------------------------------------
#Figure A5
#----------------------------------------------------------------------------------------

load('synth.Rdata')

graphics.off()
par("mar")
par(mar=c(1,1,1,1))
pdf('FigureA5.pdf', height = 8*2, width = 6*2)

  par(mfrow=c(4, 2), oma = c(0.001, 0.001, 0.001, 0.001))#
  
  title_af = paste('CAMBODIA', '\nPost-Treatment Gap (avg): ', deci(post_Gap(synth.milper.Cambodia, data.milper.Cambodia, transform=F), 5), '\n Pre-Treatment RMSPE: ', deci(pre_RMSPE(synth.milper.Cambodia, data.milper.Cambodia, transform=F), 5), sep = '')
  path.plot.cus(dataprep.res = data.milper.Cambodia, synth.res = synth.milper.Cambodia, tr.intake = 1975, Ylim = c(0, 0.07), Main = title_af, Ylab = 'Military personnel per 1,000 citizens')
  abline(v=0, lty=2)
  
  title_af = paste('CHINA', '\nPost-Treatment Gap (avg): ', deci(post_Gap(synth.milper.China, data.milper.China, transform=F), 5), '\n Pre-Treatment RMSPE: ', deci(pre_RMSPE(synth.milper.China, data.milper.China, transform=F), 5), sep = '')
  path.plot.cus(dataprep.res = data.milper.China, synth.res = synth.milper.China, tr.intake = 1949, Ylim = c(0, 0.07), Main = title_af, Ylab = 'Military personnel per 1,000 citizens')
  abline(v=0, lty=2)
  
  title_af = paste('CUBA', '\nPost-Treatment Gap (avg): ', deci(post_Gap(synth.milper.Cuba, data.milper.Cuba, transform=F), 5), '\n Pre-Treatment RMSPE: ', deci(pre_RMSPE(synth.milper.Cuba, data.milper.Cuba, transform=F), 5), sep = '')
  path.plot.cus(dataprep.res = data.milper.Cuba, synth.res = synth.milper.Cuba, tr.intake = 1959, Ylim = c(0, 0.07), Main = title_af, Ylab = 'Military personnel per 1,000 citizens')
  abline(v=0, lty=2)
  
  title_af = paste('IRAN', '\nPost-Treatment Gap (avg): ', deci(post_Gap(synth.milper.Iran, data.milper.Iran, transform=F), 5), '\n Pre-Treatment RMSPE: ', deci(pre_RMSPE(synth.milper.Iran, data.milper.Iran, transform=F), 5), sep = '')
  path.plot.cus(dataprep.res = data.milper.Iran, synth.res = synth.milper.Iran, tr.intake = 1979, Ylim = c(0, 0.07), Main = title_af, Ylab = 'Military personnel per 1,000 citizens')
  abline(v=0, lty=2)
  
  title_af = paste('NICARAGUA', '\nPost-Treatment Gap (avg): ', deci(post_Gap(synth.milper.Nicaragua, data.milper.Nicaragua, transform=F), 6), '\n Pre-Treatment RMSPE: ', deci(pre_RMSPE(synth.milper.Nicaragua, data.milper.Nicaragua, transform=F), 5), sep = '')
  title_af = paste('NICARAGUA', '\nPost-Treatment Gap (avg): ', '-0.0000279', '\n Pre-Treatment RMSPE: ', deci(pre_RMSPE(synth.milper.Nicaragua, data.milper.Nicaragua, transform=F), 5), sep = '')
  path.plot.cus(dataprep.res = data.milper.Nicaragua, synth.res = synth.milper.Nicaragua, tr.intake = 1979, Ylim = c(0, 0.07), Main = title_af, Ylab = 'Military personnel per 1,000 citizens')
  abline(v=0, lty=2)
  
  title_af = paste('RWANDA', '\nPost-Treatment Gap (avg): ', deci(post_Gap(synth.milper.Rwanda, data.milper.Rwanda, transform=F), 5), '\n Pre-Treatment RMSPE: ', deci(pre_RMSPE(synth.milper.Rwanda, data.milper.Rwanda, transform=F), 5), sep = '')
  path.plot.cus(dataprep.res = data.milper.Rwanda, synth.res = synth.milper.Rwanda, tr.intake = 1994, Ylim = c(0, 0.07), Main = title_af, Ylab = 'Military personnel per 1,000 citizens')
  abline(v=0, lty=2)
dev.off()

pst_gaps <-c(
  post_Gap(synth.milper.Cambodia, data.milper.Cambodia,3),
  post_Gap(synth.milper.China, data.milper.China),
  post_Gap(synth.milper.Cuba, data.milper.Cuba),
  post_Gap(synth.milper.Iran, data.milper.Iran),
  post_Gap(synth.milper.Nicaragua, data.milper.Nicaragua),
  post_Gap(synth.milper.Rwanda, data.milper.Rwanda)
)

